library(openxlsx)

'
描述：gee 分析 varsy[j] ~ (varsx[i] + varsc)，分别为因变量、自变量、协变量
依赖: openxlsx, geepack, gtsummary
'
gee_run = function (
  data,                # 数据框
  varsx,               # 自变量（暂时只支持变量名形式，后续会支持 index 形式） 
  varsy,               # 因变量（暂时只支持变量名形式，后续会支持 index 形式）
  varsc = NULL,        # 协变量（暂时只支持变量名形式，后续会支持 index 形式）
  file = "output.xlsx",# 输出文件名
  ...                  # 其他参数将会传递给 geeglm
) {
  # 初始化 workbook
  wb = createWorkbook()
  extract_data = data
  fit_table = list()
  for (j in 1:length(varsy)) {
    fit_table[[j]] = list()
    for (i in 1:length(varsx)) {
      fit_table[[j]][[i]] = list()
    }
  }
  # 外循环：控制因变量的个数
  for (j in 1:length(varsy)) {
    out = NULL
    ######## 统计一些数量信息
    total_obs = nrow(data)
    total_unit = length(unique(data$id))
    ########
    # extract_data = extract_data[complete.cases(extract_data[, c(varsy[j], varsc)]), ]  # 去除因变量、协变量中包含 NA 的行
    ######## 统计一些数量信息
    after_cy_obs = nrow(extract_data)
    after_cy_unit = length(unique(extract_data$id))
    ########

    # 内循环：控制自变量的个数
    for (i in 1:length(varsx)) {
      # extract_data = extract_data[complete.cases(extract_data[, varsx[i]]), ]   # 去掉自变量中包含 NA 的行
      ######## 统计一些数量信息
      after_x_obs = nrow(extract_data)
      after_x_unit = length(unique(extract_data$id))
      ########
      y = varsy[j]                                             # 拼接因变量
      x = paste(c(varsx[i], varsc), collapse = " + ")          # 拼接自变量
      formula = as.formula(paste(y, "~", x))
      print(formula)
      fit = geeglm(formula = formula,     # 使用 as.formula() 生成公式
                   data = extract_data,
                   id = id, family = "binomial",
                   corstr = "independent", scale.fix = TRUE)
      out1 = summary(fit)$coefficients
      out1$variable = varsx[i]
      out1$`Pr(>|W|)` = round(out1$`Pr(>|W|)`, 4)
      ### OR, CI_OR
      out1$OR = round(exp(coef(fit)), 2)
      CI = round(as.data.frame(exp(confint.default(fit))), 2)
      CI$OR = paste("(", CI$`2.5 %`, ", ", CI$`97.5 %`, ")", sep = "")
      out1$CI_OR = CI$OR
      ###
      ###### 写入一些备注信息
      out1$total_unit = total_unit
      out1$after_cy_unit = after_cy_unit
      out1$after_x_unit = after_x_unit
      out1$x = varsx[i]
      out1$formula = paste(y, "~", x)
      ######
      # 追加结果：内循环第一次执行时 out 为空，会直接赋值为 out1；以后都是追加行
      out = if (is.null(out)) out1 else rbind(out, out1)
      ### summary-plot
      fit_table[[j]][[i]] = fit
    }
    # 输出为单个表格
    addWorksheet(wb, sheet = varsy[j])
    writeDataTable(wb, sheet = varsy[j], x = out, rowNames = T)
  }
  # 输出到 xlsx 文件
  saveWorkbook(wb, file, overwrite = T)
  return (fit_table)
}